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We argue that the binding between doubly occupied (doublon) and empty (holon) sites governs 
the incoherent excitations and plays a key role in the Mott transition in strongly correlated Mott- 
Hubbard systems. We construct a new saddle point solution with doublon-holon binding in the 
Kotliar-Ruckenstein slave-boson functional integral formulation of the Hubbard model. On the half- 
filled honeycomb lattice, the ground state exhibits a continuous transition from the semimetal to an 
antiferromagnetic ordered Slater insulator with coherent quasiparticles at Uci — 3.4f, followed by a 
Mott transition into an electron-fractionalized AF* phase without coherent excitations at Uc2 — 5.7t. 
We show that doublon-holon binding unites the three important ideas of strong correlation: the 
coherent quasiparticles, the incoherent Hubbard bands, and the deconfined Mott insulator. 

PACS numbers: 71.10.-w, 71.27.-|-a, 71.10.Fd, 73.30.+h 



The fundamental theoretical challenge of the strong 
correlation problem is the description of both the low 
energy coherent quasiparticles (QPs) and the higher en- 
ergy incoherent excitations, as well as the spectral weight 
transfer from coherent to incoherent excitations with in- 
creasing correlation strength. Two very important ideas, 
the emergence of two broad incoherent features known 
as the Hubbard bands and the existence of rcnormal- 
ized QPs with a Luttinger Fermi surface (FS) were ad- 
vanced by Hubbard [l|, and Brinkman and Rice re- 
spectively. Unfortunately, the Hubbard equation of mo- 
tion scheme that produces the incoherent spectrum fails 
to produce QPs correctly and violates Luttinger's theo- 
rem ; whereas the approaches based on the Brinkman- 
Ricc-Gutzwiller wave functions Q capture a Luttinger 
FS of QPs, but find serious difficulties in constructing 
variational excited states to account for the incoher- 
ent excitations. Faced with this enigma, numerical ap- 
proaches such as exact diagonalization, quantum Monte 
Carlo (QMC), and the dynamical mean field theory 
have played a key role in recent years. In this paper, we 
develop new analytical insights, construct a unified the- 
ory for both coherent and incoherent excitations and the 
Mott transition, and study the Hubbard model on the 
honeycomb lattice in view of the recent debate over the 
emergence of a gapped spin hquid (SL) phase 

For simplicity, we consider the half-filled single-band 
Hubbard model. The local Hilbert space consists of 
empty (holon), doubly occupied (doublon), and singly 
occupied sites. The Brinkman-Ricc-Gutzwillcr approach 
amounts to a metallic state where the holon, denoted 
as a boson e, and the doublon, as d, condense fully 
with macroscopic phase coherence, as obtained by the 
Gutzwiller approximation ll| or the slave boson mean- 
field theory (SBMFT) [Uy^The metal-insulator transi- 



tion is thus forced to follow a route where the density 
of doublons and holons vanish together with the conden- 
sates: rid = rie = (d) = (e) = such that there is exactly 
one electron per site. As a result, single-particle motion, 
coherent or incoherent, is completely prohibited. This 
so-called Brinkman- Rice (BR) transition is different from 
the Mott transition induced by the complete transfer of 
the coherent QP weight into the incoherent background, 
i.e. the depletion of the condensate while keeping the 
doublon/holon (D/H) density an analytic function in the 
correlation strength U. 

We will show that the crucial physics uniting the dis- 
parate ideas of Hubbard and Brinkman- Rice is the bind- 
ing between doublon and holon: {diCj) ^ 0. In the Mott 
insulator at large U, although the D/H condensate van- 
ishes ((d) = (e) = 0) together with coherent QP, the 
D/H density remains nonzero (rid = n-e ^ 0). The in- 
coherent motion of the QP is thus possible by breaking 
the doublon-holon (D-H) pairs, giving rise to the higher- 
energy incoherent excitations. With decreasing U, the 
D/H density increases and the D-H binding energy de- 
creases. At a critical Uc, the D-H excitation gap closes 
and the D/H single-particle condensate develops, mark- 
ing the Mott transition. On the metallic side of the Mott 
transition, D-H binding continues to play an important 
role since an added electron can propagate either as a 
coherent QP via the D/H condensate or incoherently via 
the unbinding of the D-H pairs. 

The idea that D-H binding plays a role in Mott- 
Hubbard systems was introduced by Kaplan, Horch, and 
Fulde [l3| and studied in the context of improved vari- 
ational Gutzwiller wave functions [13, [3l ■ It has been 
difficult to advance because of the difficulty in the calcu- 
lations and particularly in constructing the appropriate 
wave functions for excitations. We will show that the 
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FIG. 1: Schematic phase diagram obtained for the half- filled 
honeycomb lattice Hubbard model with D-H binding. 

physical picture discussed above can be realized in the 
slave-boson functional integral formulation of the Hub- 
bard model introduced by Kotliar and Ruckenstein (KR) 
12| by constructing new saddle point solutions that in- 
clude the D-H binding. This approach also offers a treat- 
ment of the magnetism at half-filling that compares well 
with QMC simulations [iBl , and has the added advantage 
of allowing the study of excitations and finite tempera- 
ture properties fl7j . 

As a concrete example, we studied the D-H binding 
in the half-filled Hubbard model on the honeycomb lat- 
tice at zero temperature and obtained the phase diagram 
shown schematically in Fig. 1. A continuous transition 
from the semimetal (SM) to an antiferromagnetic (AF) 
ordered insulator takes place at a critical Ud — 3.4<, 
suggesting that the gapped SL phase proposed by Meng 
et. al. 16| may correspond to an AF ordered phase in 
the thermodynamic limit. Sorella et. al. [lo| recently 
extended the QMC and the finite size scaling analysis to 
much larger system sizes and discovered that the signa- 
ture of the gapped SL disappears and is replaced by that 
of a continuous SM to AF order transition at J7 ~ 3.8t, 
in qualitative agreement with our results. Remarkably, 
we found a second quantum phase transition at a critical 
Uc2 — 5.7i beyond which the D/H single-particle conden- 
sate vanishes (do = 0) amid a finite density of doublons 
bound to holons. For U > Uc2, a new AF phase, termed 
as AF* in Fig. 1, emerges where the electrons are frac- 
tionalized and there are no coherent QP excitations. 

In KR's theory, the electron operator is written as 



Ltaielpia- + Pigdi)Riaficy, whcrc the boson oper- 
ators describe the holon (e^), doublon (di), and singly- 
occupied sites (pia), and fi^- is a fcrmion operator. The 
operators and i?^ are diagonal with unit eigenval- 
ues in the (empty, a) and the (cr, doubly-occupied) sub- 



spaces, respectively, Lia- = {l—d]d. 



■Y 



and Ria = 

(1 - e[ei - pl^p^sY■ KR found that for a = ;3 = -1/2, 
the saddle point solution recovers the Gutzwiller results. 
The Hubbard model Hamiltonian is thus given by [l2[ , 



where Zi„ = Lia{e\pi„ + pYdiJRia and the Lagrange 



multipliers \i and Xia enforces the local constraints: 
Y.aPlcrPi'y + d\di -1 = and Qia = 
di — 0. The partition function is a 
functional integral over the quantum fields 22|: Z = 
J !?[/, e,p, d. A] cxp(— Jjj^ £fir), whcrc the Lagrangian 

^ = Eii4dreMdrd,) + E.ApldrP^. + fldr.M + H. 

The KR saddle point corresponds to condensing all bo- 
son fields uniformly in Eq. ([T|). The results on the hon- 
eycomb lattice [23| can be summarized: Restricting to 
the PM phase, the doublon density dg decreases linearly 
from 1/4 at [/ = and vanishes at the BR metal-insulator 
transition Ubr — 12.6t. When magnetism is allowed, a 
SM to an AF insulator transition arises at Um — 3. It. 
The D/H condensate remains nonzero for any finite U 
and the AF phase is a coherent Slater insulator 24| . 

To construct the new saddle point solution with D- 
H binding, we need to include the dynamics of the 
bosons. Introducing the operators for the D-H pairing 
Aij = diCj, and the D/H hopping xfj = d\dj, x% = e^j 
on the nearest-neighbor bond, as well as the density 
operators fif = djd;, hf = eje^ on each site, we can 



write in Eq. H]): zja^ja = 9tj[{x%Yp\aPjo 
^y^aPJal whcrc .g^ - {Yf.V 



^TiP\aP]a 



1/2 Y° 



( 1 -plsPif ) ( 1 - PjaPiT ) - n1 ( 1 ~p]aPj<y ) - ( 1 -plsPif ) + 
|Ajip. These inter-site correlations can be treated explic- 
itly by writing the partition function as a path integral 
over the corresponding fields and taking the saddle-point 
solution [2J|. The results can be understood intuitively 
as minimizing the variational Hamiltonian obtained by 
substituting the above expression for zj^zj^ into Eq. ([1]) 
with the operators replaced by their expectation values, 
and adding the corresponding variational terms: 



H 



DH 



- t E + h.c.) + uY,dld, 

i i,(T ia 

-iT.\^f(4d,-4)+x:f{4^,-x^,) 

(id) 

+ Kjid^ej - A,j) + AJ,(e,dj - Aj.) + h.c. 
+ z5][6f-(dld.-nf) + er(ele,-<)], (2) 

i 

d,v e,v d,v 



where A^^-, x°/, xl'j' , <^T^ , and e^'" are self-consistently 
determined variational parameters. Note that the pia- 
bosons are treated as condensed fields for simplicity since 
their density (single occupation) is large (25j . 

We next discuss the saddle point solution of Eq. 
on the honeycomb lattice with 2N sites. Symmetry re- 
quires the expectation values Aij = Ad, xfj = xtj = Xd, 



^Xi 



Ud, and correspondingly, iAJj = AJJ, iXif 
'^'^ = Cj. Moreover, iXi = A, iXao 



V ■ f^^v 

Xd, ^^i 
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AY) 



PbI- — J^Ocr, where A and B denote 



i^Bg = Act, and p)l 
the two sublattices. The normahzation factor gfj = g = 

withF, = l-2na-2pl, + 2pl,nd+pt + Al 
As shown later, this expression recovers the noninter act- 
ing hmit at U = 0. The Haniihonian thus simphfies to 
^DH = Hf+Hd + l2Nixlxd + A^Arf) - 2iV(e^ + 2A) - 
ANe^Ud + '2NJ2ai^ ~ ^(T)Pocr- The fermion spectrum is 
given by 



k,o- 





t - 







fAka 
fskcr 



, (3) 



where ?7k = exp(ifcj,) + 2 cos(V3fci:/2) exp(— ifcj,/2), ~ 

tg 2potPoiXd + (Pot + Pli)^d , and the chemical poten- 
tial /i = C//2 by the particle-hole symmetry. The fermion 
dispersion is thus £'±(k) = ±^{\^ - A;)2/4 + |x/?7kp- 

A fermion gap = jA-]- — Aj^| would open in the Dirac 
spectrum with AF order. Similarly, the charged boson 
spectrum is determined by Hd = X^k ^k-^^k^k, where 
^'k = [dAk,rfBk,e^j^,e^j^]^ and 



£d 
"Xd'/k 





-XdVk 

^d 







-XlVk 





-XdVl 
el 



(4) 



Here 



A, el ^ ~ig\fX)T.A^-ph)Y.. and 

Xl = 2tgpotPoiXf, K = 9XfY.a{tpla - 9^dX}Y<^)- 
The boson dispersion is obtained by Bogoliubov trans- 



formation: (k) 
the condition > 3(|Xd 



|A-7?kp under 
The boson energy 



gap is thus Ed - 2^(e- - 3\xl\r - (3|A-|)^ It is in- 
structive to examine the equations for the D/H density, 
D/H hopping, and the D-H binding. 



nd=dl + ^ ^idl^^dc^k) 

a={A.B} k 



(5) 



iXd, Ad) =dl + ^ J2^Vkd\^{dBk, e^k) + ^-C-)- (6) 

k 

The closing of the boson gap Ed leads to a zero energy 
mode at k = whose occupation enables the single-boson 
condensate dp, and is subsequently taken out of the mo- 
mentum summations. 

The ground state is determined by solving the self- 
consistency equations derived from minimizing the en- 
ergy IH]: {dHHo/dxi) = and x = {Ad, Xd, nd, poa, 
Act, A). Once the saddle point solution is obtained, the 
spectral function can be obtained from the one-electron 
Green's function Gccr(k,T) = -{TrCaka{T)cli^^{0)). 
Since the latter involves convolutions of the d/e boson 
normal and the anomalous (due to pairing) Green's func- 
tions with those of the /^--fermion, the single-particle 



energy gap for the physical electron is the sum of the 
fermion and boson gaps E = Ed + Ef. More importantly, 
the coherent QP excitations would only emerge with the 
D/H condensate that recombines the charge and spin de- 
grees of freedom, which can be detected by the coher- 
ent peaks in the integrated spectral function (ISF) 
Na{uj) = — Im dre^'^^ X^k a G'ao-(k, t). Our results are 
plotted in Figs. 2 and 3. The parameter A^ increases 
with U and the D-H pairing field A^ is nonzero for all 
finite-[/ (Fig. 2a). At f7 = 0, all doublons and holons are 



single-particle condensed dp 



pIo 



1/4, recov- 



ering the noninteracting limit. The SM phase remains 
stable at small U . With increasing [/, the doublon den- 
sity decreases. Due to the increase in D-H binding, the 
D/H condensate do decays faster than in the SBMFT. 
Let's first restrict the solution to the PM phase with 
Pot = Poi- As shown in Fig. 2b, there is a Mott tran- 
sition at Uc = 7.3t, considerably smaller than 12.6t for 
the BR transition. The condensate do vanishes and all 
doublons are bound with the holons in the Mott insu- 
lating phase for U > Uc, accompanied by the opening 
of a charge gap Ed linear in [/ — t/c(Fig. 2c). The ISF 
is shown in Fig. 3a. Notice the transfer of the coherent 
QP weight to the incoherent part with increasing U and 
the complete suppression of the coherent QPs in favor of 
two broad incoherent spectra beyond the Mott transition 
from the two branches of bosonic excitations (k) [^J] . 
Since the /icr-fermion spinon remains gapless, the insu- 
lating phase is a gapless SL. We do not see evidence for 
the proposed gapped SL phase [1] . 

Next, we allow magnetism and study the interplay be- 
tween AF order, D-H binding, and the Mott transition 
in the ground state. Fig. 2b shows the onset of staggered 
magnetization (m) accompanying the AF transition at a 
critical Ud — 3.4t, which compares well to the contin- 
uous SM to AF transition observed in the most recent 
QMC calculations on large systems [l^l- We find that 
for Uci < U < Uc2, where Uc2 — 5.7t, although a single- 
particle gap Ef opens in the fermion sectors (Fig. 2c), 
the zero energy mode remains stable in the d-e sector and 
continues to support the D/H condensate. Thus, the spin 
and charge recombine in this regime and there are coher- 
ent excitations corresponding to the sharp QP peaks in 
the ISF shown in Fig. 3b. This phase is thus an AF Slater 
insulator. Remarkably, a Mott transition takes place at 
Uc2- For U > Uc2, a new AF phase, the AF* phase, 
emerges with the opening of the gap Ed ccU — Uc2 in the 
d-e sector (Fig. 2c), where the D/H condensate vanishes 
as all doublons are bound to holons. Thus electrons are 
fractionalized in the AF* phase and there are no coher- 
ent excitations, as shown in the broad ISF in Fig. 3b at 
U = 6t. The gap for the physical electron, E = Ef + Ed, 
exhibits a derivative discontinuity at Uc2 (Fig. 2c) due to 
the opening of the Mott gap Ed in the AF* phase. Note 
that the magnetization m is analytic across the AF^-AF* 
transition (Fig. 2b), which is a topological confinement- 
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FIG. 2: (color online) Plotted as a function of Hubbard TJ 
in (a) are variational parameter for d-e binding and D-H 
pairing order parameter A^; in (b) D/H density n^, conden- 
sate density dgi s-nd AF staggered magnetization m; and in 
(c) ground state energy per site E and energy gaps in the 
boson and fermion H/ sectors. The corresponding results 
in the restricted PM phase are also shown. 

deconfinement transition associated with the Ising-like 
global Zi symmetry (di — > — di, — > — e^) that is broken 
in the AF phase by the D/H condensate and restored in 
the AF* phase. Finally, let's examine the saddle point 
stability with respect to gauge field fluctuations. It is 
known that KR formulation introduces three J7(l) gauge 
fields [26j since the action is invariant under: — >■ 6^6*^' , 

and Xi + Oi, X^a Xia + - 4>ia- The p^^ con- 

densate breaks two of the U(l) symmetries and turns the 
gauge fields associated with (j)icr massive by Anderson- 
Higgs mechanism. The remaining U{1) symmetry is also 
broken in the SM and the AF phase by the D/H con- 
densate, making the ^i-gauge field massive. In the AF* 
phase, it is the D-H pairing Ay that breaks the U{1) 
symmetry and the ^^-gaugc field remains massive, as does 
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FIG. 3: (color online) Integrated spectral function showing 
the transfer of coherent QP weight to the incoherent spectra 
with increasing U listed in the legends, (a) PM phase: Mott 
transition at Uc — 7.3t. (b) Ground state: SM to AF transi- 
tion at Uci ~ 3At and AF to AF* transition at Uc2 ^ 5.7t. 



its staggered component due to the D/H hopping fields 
Xif- Tlie absence of gapless U(l) gauge field fluctuations 
supports the stability of the obtained phases. 

In summary, we have shown that the D-H binding plays 
an essential role in describing the incoherent excitations 
and the Mott transition in strongly correlated systems. 
For the honeycomb lattice Hubbard model, wc showed 
that the SM to AF Slater insulator transition is followed 
by a Mott transition into a fractionalized AF* phase with 
increasing U. Interestingly, a different AF* phase of a 
fractionalized antiferromagnet was proposed in the effec- 
tive Z2 gauge theory description of doped Mott insulators 
in the projected Hilbert space (U = 00) where spinous 
are paired into a Neel state and doublons are absent 27 1 . 
In contrast, the incoherent charge excitations through D- 
H binding is essential in the AF* phase proposed here, 
which is more inline with the importance of doublons in 
describing Mottness emphasized recently [28|. Such an 
AF* phase on the square lattice may have been observed 



in the parent compound of the high-Tc cuprates by angle- 
resolved photoemission [2^ . 
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99ER45747 and de-sc0002554, and the 1000 Talents Plan 
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SUPPLEMENTARY MATERIAL 



Following the discussions in the main text, the partition function for the Hubbard model can be written as a 
functional integral over the quantum fields of the spin-carrying fermions fi^ and the slave bosons e;, d;, and pia 
introduced by Kotliar and Ruckenstein (KR) P , 

Z=J V[f, p]V[e,e^V[p,p^V[d,d^]V[X,X„]e- ^0 "^"^ (SI) 



' where H is given in Eq. (1), 



The functional integral over the Lagrange multipliers and A^o- enforces the following local constraints 

= e\e, + p\aP^<y + d\d,~l = 0, (S4) 



V ■ /I/- - pIp- - - 0, Va. (S5) 

' The electron operator is written as Ci„ — Zi„Jia where Zi„ ~ Li„{e\pia -\- p\^di)Ria with the normalization factors 
^ . Ria = (1 — ejei ~ pl^pis)"^/"^ and Li„ = (1 — d\di ~ pl^pia)^^^^ ■ Because of the projection-operator property of the 
— ' ' bosons, the normal ordering of the square roots in the hopping renormalization factor zj^Zia- is not an issue, as they 
can be translated directly into c-number boson fields in defining the functional integral. 

Q ' A. Conventional slave-boson mean-field theory on honeycomb lattice 

o . 

In the saddle-point approximation of the action considered by KR [l| , which is referred to as the conventional slave- 
boson mean-field theory (SBMFT) or the Gutzwiller approximation, all the boson fields in Eq. (jSSp are condensed: 
' Ci ~ eo, di ~ do, and pia = Poa- For the honeycomb lattice Hubbard model at half- filling, the SBMFT was first 
, studied by Fresard and Doll ■ For the purpose of comparison, we plot the solutions in Fig. [STl which should be 
■ contrasted to the new saddle point solution in the presence of doublon-holon (D-H) binding shown in Fig. 2 in the 
l/^ , main text. In the PM phase, the doublon density dp decreases linearly from 1/4 at [/ = and vanishes at Ubk — 12. 6t, 
' corresponding to the Brinkman-Rice metal- insulator transition fFig. lSlb ). When magnetism is allowed, the semimetal 
, to AF insulator transition takes place at Um — 3. It with the onset of the staggered magnetization m = p^^ — p^^. 
I The doublon/holon (D/H) condensate density cZq = e§ remains nonzero for any finite Hubbard U, so the AF ordered 
phase in the SBMFT is a Slater insulator with coherent QP excitations above the AF insulating gap. The latter is 
determined by the energy gap in the fermion sector 'E.f and scales linearly with U — Um as shown in Fig. ISlb . 

, B. Doublon-holon binding and intersite correlations of the slave bosons 

In order to construct the new saddle-point solution that takes into account the effects of D-H binding, we need to 
consider the dynamics of the bosons and the their intersite correlations. To do so, we introduce the operators for D-H 
pairing Aij = diCj and the D/H hopping xfj = dldj, Xij = ejej on the nearest- neighbor bonds, as well as the H/D 
density fif — d|c?i, hf — e\ei on each site. The hopping renormalization factor in Eq. (|S3p can thus be written as 



(Xy- ) VLPj> + xfjPiap}^ + ^^i^p\aP]s + ~^\jV^sP3a\ , (S6) 

where the prefactor = R}^^L\^LjaRj(T can be decoupled according to 

9^ = ^ = (1 - - v]aVio) - nl(\ - pIp,.) ~ n^{l - + | A,, \\ (S7) 
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FIG. SI: (color online) Conventional SBMFT solution of the Hubbard model on half-filled honeycomb lattice. The Hubbard 
U dependence of (a) condensate density d^, staggered magnetization m, and (b) energy per site E and single-particle energy 
gap in the fermion sector 5/. 



Note that the above expression for gf^ involves exphcitly the D-H pairing, but not the H/D hopping operators. As 
we shall show later, this choice of the decoupling ensures the correct description of the noninteracting limit at [/ = 
where {Lij) = dl ^ 1/4 and (x^^) = = (Xy) = do such that z^^z^^ = 1. 

The obvious challenge is how to build these correlations into the calculation of the partition function. Since they 
enter through the rather formidable factor zj^Zj^, the usual procedure of introducing the corresponding correlation 
fields (Ay, X 



Xij, nf, nf) as Hubbard-Stratonovich fields in the path-integral docs not work here. Wc found 
that the difficulty can be overcome by introducing in the functional integral additional Lagrange multipliers in the 



corresponding channel, AJ" , Xi 



and 



such that the partition function becomes: 



DH, 



(S8) 
(S9) 



where 



.^[ef^(dld,-nf) + er(ele,-<y 



(SIO) 



The effective H-D binding Hamiltonian in Eq. (|S10p is just the one given in Eq. (2) of the main text, where the factor 
ZiaZja has the operator form given in Eqs (|S6[) and (jS7p but with the bosonic operators replaced by the corresponding 
correlation fields (A^ , xf^ , xfj j "-f i A few remarks are in order. (1) Eqs ((S8)) - (|S10|) provide an exact representation 



of the Hubbard model; carrying out formally the last two functional integrals in Ec^. (|S8p recovers the KR formulation 
given in Eqs (|S1I) - (|S3[ ). (2) In the above derivation, the intersite correlations of the bosons were not included 
explicitly for simplicity. The reason is that the pio- bosons represent single-occupation; their density is large at half- 
filling. Thus, the "pia bosons can be treated appropriately as condensed fields. We have verified that including their 
intersite correlations only leads to small quantitative changes in the results. (3) From the perspective of finding a 
saddle point solution of the action, the effective Hamiltonian iJnD in Eq. (|S10p can be understood intuitively as a 
variational Hamiltonian describing the effects of intersite correlations of the doublons and holons, including that of 
H-D binding, as stated in the main text. 
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C. D-H binding saddle point solution and self-consistency equations 



The saddle point solution corresponds to approximating the path integral in Eq. (|S8[) with a particular configuration 
of the quantum fields that minimizes the action. We consider here translation invariant paramagnetic (PM) and two- 
sublattice antiferromagnetic (AF) saddle point solutions on the honeycomb lattice at half-filling. The bond variables 
are real and isotropic and the symmetry requires A^^^ = A^, xfij) = X{ij) = Xd, = nl = Hd, and correspondingly, 

lA^.^.^ = A^, ix^-J) = «X(iJ) = Xd' " = *<^r'' = <^d- Moreover iXi = A, iXacf = i^Ba = A^-, and pAOa- = Pbos = Poa, 
where A and B denote the two sublattices on the honeycomb lattice. The normalization factor in Eq. (|S7p takes on 
the following saddle point value on the nearest neighbor bonds, 



5fe> = .9 = n ^itli = 1 - 2nd - 2pl^ + 2pl„nd + pt + 

a 

Substituting these quantities into Eq. (jSlOp . we obtain the saddle point Haniiltonian in the main text, 
H'^^ ^Hf + Ha + l2N{xlxd + A^A^) - 2N{el + 2A) - 4iVe>d + 2iV^(A - K)pI„ 



(Sll) 



(S12) 



on the honeycomb lattice with 2N number of sites, where Hf and Hd describe the energy spectrum in the fermion 
and boson sectors. In momentum space, the fermion part is given by 



fAka 


t - 







where the A-B sublatticc hopping amplitude xj 



-X/'?k - l-i 



'^PotPoiXd + ipl^+pl{)Ad 



and r]k 



2 cos(-\/3fca;/2) exp(— ifcy/2) on the honeycomb lattice. At half-filling, the particlc-hole symmetry requires 

= and At + A| = U. 



(S13) 
(S14) 



We can thus write Ao- = /i — ere, where e = {X^ — Xi)/2 is equal to zero in the PM phase, but nonzero in the AF phase. 
The fermion action can therefore be diagonalized by the transformation 



(S15) 







"Ako- 


VAkcr 




la- 


_(k)e 


~T£;{(k) ' 






"Bko- 








-(k)e 





where 7CT±(k) are the fermion quasiparticles (QPs) with eigenstate energy dispersions: 



lx}^k 



The corresponding eigenvectors, satisfying UAv.a ~ ^'sko- ~ "^^ko- and UAko- = ^""sko- ~ ^kcr, are determined by 



2 1 



El 



Wko 



1 1 



ere 



E. 



MkcrWkcr 



Xf^ 

2EI ■ 



(S16) 



(S17) 



The fermion spectrum in Eq. (|S16p shows that in the semimetallic PM phase (e = 0), the QP has Dirac dispersion 
with a renormalized velocity xj; whereas a fermion energy gap 



S/=2|e| 



(S18) 



opens at the Dirac points in the AF ordered phase. The fermion hopping and density per spin are readily obtained, 

X/ =X/ 



~ ^(/iker/^kCT) = ]y E^-^Bkff/skff) = Tp^ E \^ ) 
k k k \ -^k / 



(S19) 
(S20) 



4 



The spectrum of the charged (rf-e) boson sector is determined by Hd in the saddle point Hamiltonian (jS12[) . which 
can be written in a 4 x 4 matrix form in momentum space, 



d-Ak 
dsk 

J 







-xWk ^d + A -AX 
-Kv*k el + X -xWk 
-A;j77k -xliVk + \ 





dAk 




dsk 




- 








p - 



(S21) 



where the relations due to particle-hole symmetry in Eq. ()S14[) have been applied. The bosonic action can be 
diagonalized by the following Bogoliubov transformation, 



" dAkir) ' 
dBkir) 




UAk+ 


UAk- 


VAk+ 


VAk 




UBk+ 


UBk- 


VBk+ 


VBk 




VBk+ 


VBk- 


UBk+ 


UBk 


1 4(-) \ 




_ VAk+ 


VAk- 


UAk+ 


UAk 



a+(k)e--^+('^) 

-T_El{k) 
rE^ik) 



a_(k)e 



6T_(k)e- 



(S22) 



where a± (k) and b± (k) are four bosonic QPs with two branches of doubly degenerate eigenvalues 



(S23) 



under the condition + A > 3(1x^1 + I ^dl) ^^^^ ensures the eigenvalues to be real and positive for bosons. The factor 
of 3 in the latter arises from the maximum value of 77k- The energy gap in the d-e boson sector is thus given by 



Srf = 2 



V/(^S + A- 31x^1)2 -(3|Av|)2 



(S24) 



Analytical expressions of the corresponding eigenvectors {u, u} are more lengthy than informative and thus omitted 
here. In practice, they are computed numerically. 

The saddle point solution for a given Hubbard U can be determined by solving the set of self-consistency equations 
derived from minimizing the energy with respective to the ten variables {xd. Ad, rid, poo-, £, A, Xdi ^d }• 



{dHuF/dxd) 
{dHuF/dAd) 



Xd = 2i,9mmx/> 

A2 = gX/E(^Po.-3AdX}>^- 



{dHMv/dud) = 

{dH,s/dpna) = 
(dHee/de) = 
(dHes/dX) = 

{dHMF/dxD = 
(9ffMF/5A2) = 
(dHuF/de^i) = 



ed = -3,9'x/X/E(l-Po*)^- 

(J 

[A - U/2 + ere - 6g^x/X/^ff (1 - nd - Pq^)] poa = QtgxfixdPos + AdPoa), 

2nd+Pot+Poi = 1' 

Xd = dl + {Vkd^i.dBk + h.c.), 

k 

Ad = rfo + ^ X! ' (^kdAkSBk + h-c.) , 



nd = dl + ^ X X ' {dakdck) 

a = {A,B} k 



(S25) 

(S26) 

(S27) 

(S28) 
(S29) 
(S30) 

(S31) 
(S32) 
(S33) 



It is instructive to examine the last three equations for the D/H hopping, D-H pairing, and the D/H density fields. 
The closing of the gap Sd in the boson spectrum gives rise to a zero energy mode at the F ~ (0, 0) point [see Fig. ((S3l) ] 
whose occupation enables the single-boson condensates do ~ ^ 0. Thus, whenever do 7^ 0, the F point is removed 
from the momentum summations in Eqs (jS311 IS321 IS33P such that the remaining summation over k comes solely from 
the contributions introduced by the D-H binding. Accordingly, the solutions to this set of self-consistency equations 
(|S25IIS33l) must be searched with two attempts: (i) assume o?o = and solve these ten equations for the ten unknown 
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variables and (ii) assume a nonzero do. The latter introduces one more unknown variable (do), while the required 
condition for the existence of the zero energy mode gives us the extra equation 

A = 3(1x^1 + IA2I), (S34) 

and thus one can determine the state by solving the eleven self-consistency equations for the eleven unknown variables. 
If several solutions are found for a given U, the one with the lowest energy should be chosen as the ground state. In 
practice, we found only one solution for any given U. 

Since most of the results have been discussed in the main text, we provide here only a few additional details. In 
Fig. IS21 the self-consistently determined variational parameter Xd ^^'^ the D/H hopping parameter Xd are plotted 
as a function of U in the ground state, as well as in the restricted PM phase. At U = 0, x^d 0. However, since 
= (Fig. 2a in the main text) at J7 = leads to = 1/4, the noninteracting limit is correctly recovered. The 
variational parameter Xd increases and the D/H hopping field Xd decreases monotonically with increasing U. Note 
that Xd remains finite for any finite U and evolves analytically across the Mott transitions at Uc = 7.3t in the PM 
phase and at the AF— >-AF* transition at Uc2 ~ 5.7t in the ground state, as does the D-H pairing field shown in 
Fig. 2a in the main text. As discussed in the main text, the nonzero values of Xd and ensure that both the uniform 
and the staggered components of the charge U(l) symmetry are broken and the associated gauge fields are massive 
in the electron- fractionalized AF* phase. 

It is also enlightening to discuss the boson energy spectrum and the integrated spectral function (ISF) in the 
doublon/holon sector, which are plotted in Fig. [S3] in the PM phase at several values of U across the Mott transition. 
Their behaviors in the AF and AF* phases are similar. At a given U, the boson spectrum shows the two doubly- 
degenerate dispersive branches obtained in Eq. (jS23P and Eq. (jS24p . There are several noteworthy features, (i) 
Both dispersive branches are flat around the M point of the hexagonal Brillouin zone, leading to the two Van Hove 
singularities in the boson ISF shown in the right panel of Fig. IS3I (ii) The two branches cross and produce the 
Dirac-cone-like behavior at the K point at a finite energy that increases with increasing U, leading to the V-shaped 
density of states. Remarkably, (i) and (ii) combine to form a Dirac-cone like dispersion that is similar and can be 
regarded as a "ghost" band of the bare electron dispersion carried by the excitations of the doublon/holon complex. 
This property was pointed out in the systematic large-N expansion study of the t-J model for doped Mott insulators 
Q . It is remarkable that the ghost Dirac-cone feature manifests itself in the broad peak-dip-peak structure in the ISF 
of the physical electrons shown in Fig. 3 of the main text, which can be identified as the holon-doublon excitations, 
(iii) The low energy properties of the boson dispersion near the F-point is intriguing and important. For U < Uc, i.e. 
on the metallic side of the Mott transition, the lower energy branch is gapless, i.e. ^^^(k) = and disperses linearly 
away from the F point. The existence of the zero-energy mode, together with the vanishing of the DOS/ISF N{E), 
enables the finite-temperature D /H condensation in a two-dimensional system such that 7^ at zero temperature. 
In contrast, for U > Uc, ot U > Uc2 when magnetism is allowed, a charge gap 7^ opens up at the F point. 
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FIG. S3; (color online) The energy spectrum along high symmetry directions (left panel) and the ISF (right panel) of the 
doublons and holons in the PM phases at several Hubbard U across the Mott transition at Uc = 7.3t. The contribution from 
the condensate on the semimetallic side of the Mott transition is not shown. 



indicating the emergence of the Mott insulating or the AF* phase with complete suppression of the D/H condensate. 
Note also that the gapped i?l(k) is parabolic near F, giving rise to a finite ISF N{E) at the band bottom. 

D. Integrated spectral function for the physical electrons 

We provide here additional discussions on the calculation of the integrated spectral function (ISF), which is also 
the tunneling density of states (DOS), for the physical electrons. The electron spectral function is defined in terms of 
the retarded single-particle Green's function Go-(k, w) 0]: 



Aa{k,LL!) = -Im^ Gc„^(k, w), Gao-(k, r) = -{TrCaUcr{T)cl^^^{0)). 



(S35) 



Because of the presence of two sublattices on the honeycomb lattice, the Green's function and the spectral function 
acquire an additional sublattice index a in the above equations. The ISF or the tunneling DOS, is given by 



N^iio) = - V ImG„,(k, c^) = - V Im / dre'"^G„,(k, r 
I, _ 1, _ Jo 



(S36) 



As discussed in the main text, the electron operator is composed of = Lic{e\pi„ + p\gdi)Ri„ fi„ in the KR slave 
boson formulation [J] . Within our saddle point solution, we keep the product of the fermion and slave-boson operators 
and approximate the normalization factors by the saddle point average. The electron operator in momentum space 
is given by. 



ikiT — Taa ^ ^ (^ctqPQq'iT + P^jq/g-dq) /a,k— q— q',(T j 

q,q' 



with the normalization factor 



(S37) 



(S38) 
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The electron Green's function in Eq. (|S35|) is therefore given by 

G„,(k,T) =rL5lA„.(q,q',T)GL(k-q-q',r), (S39) 

q.q' 

where is the /cr-fermion Green's function and A involves the normal and anomalous (due to pairing) Green's 
functions of the bosons 

Aoa(q,q',T) =(r,eiq(T)e„q(0))(T,p«q-,(T)p|^q,,(0)) + (T,d„q(r)dlq(0)) (r,p^,^(TKq,^(0)) 

+ (T,eiq(r)4q(0))(r.p„q.,(rKq,^(0)) + (T.d„q(r)e„q(0)) (r,p^,^(T)piq,,(0)). (S40) 
The ISF of the physical electrons in Eq. (|S36p becomes 

N^iuj) = -J2 '-Llm / dre"^"A„,(r)GL(k, r), where A„,(t) = ^ A„,(q, q', r). (S41) 

In frequency space, the spectral function thus involves convolutions of the boson Green's functions with that of the 
/cr-fermions, and as a result, the single-particle energy gap in the DOS of the physical electron is the sum of those in 
the fermion and boson sectors: S = S/ + + Sp. Within our theory, the po-boson condensate exists over the entire 
phase diagram and thus Sp = 0. 

It is instructive to write each boson operator as the sum of the condensate and fluctuations: 6^'' = ^o'^k + ^k \ where 
b stands for the {d,e,p„) bosons. Although this is not necessary, doing so facilitates well the following discussions of 
the coherent and incoherent contributions to the electron spectral function. Note that the fluctuations b'^'^ are boson 
operators, obeying boson commutation relations and the energy spectrum discussed above. Thus, the normal and 
anomalous boson Green's functions can be written as 

{Trb^^\r)b'J^^HO)) = bob'^S^S^ + {T,b<^\r)b'(^HO)). (S42) 

Decomposing the condensate and fluctuation contributions in Eq. (|S40p this way and keeping the leading order 
fluctuations involving a single boson Green's function, the Aaa{T) in Eq. (|S4ip becomes 

A„.(r)=Ar'*(r)+Ar(r), 
where the condensate part A^™'^(t) — (iQ(pot +Po4.)^i and the fluctuation part 

^aTir) = [iTrPa^a{r)pi^AO)) + (T.^q^ (r)p„q^ (0))] +pL.E(^-^~UMeaq(0)> +pLE(^-'^~"qM^~"q(0)) 

q q 

. E [(r.elq(r)dt q(0)) + (T.d„q(r)g„q(0))] . (S43) 



q 

- PotPoi 



Correspondingly, the ISF in Eq. (jS41[) can be written as 

No,{u;) = 7V-'^(w) + ^r°'X^), (S44) 

with 

^coh(incoh)(^) = -E^-Inr /^dre--Ar''(«-)(r)GL(k,r). (S45) 

The coherent part of the ISF comes from the single-boson condensates that recombine the charge and spin degrees of 
freedom, leading to coherent quasiparticle (QP) excitations associated with the coherence peaks in the ISF. 
The /o.-fermion Green's function is easy to compute in terms of the fermion QP wave functions in Eq. (jSlSp . 

GL(k,r) = -{T,.UAr)fLAO)) = ~\uc.u.\'nF{El)e^^^ - | (S46) 

where np{E) = l/[exp(/3_E) + 1] is the Fermi distribution function and j3 ~ l/ksT with T the temperature. As a 
result, the coherent contribution to the ISF is given by 

A^r'Xw) = TTdlipot+Poif E [l""k.P<5(c. + Ei) + I wP^K - e{^)] . (S47) 

k,(T 
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FIG. S4: (color online) The coherent (shaded red areas), incoherent (blue solid lines), and the total (black solid lines) integrated 
spectral function at different values of Hubbard U in (a) the restricted PM phase where the Mott transition is at Uc — 7-St and 
(b) the ground state where the AF to AF* transition is at Uc2 — 5.7t. 



Because lUakcrP + l^alcTl 



1, the total coherent spectral weight of the ISF is 



In the noninteracting limit at [/ = 0, the saddle point solution has do = poa = 1/2 and raa = 2, such that the 
total coherent weight W^°^ — 1 with the complete suppression of the incoherent part of the spectral function. In the 
opposite limit on the Mott insulator side, where the doublon/holon condensate do = eo = 0, the coherent spectral 
function is completely suppressed W™'' = and the entire spectral weight must come from the incoherent part of the 
spectral function. 

In the incoherent part of the ISF, defined in Eqs (jS45P and (jS43p . the convolution of the boson and fermion Green's 
functions gives broad spectral features. Since in our saddle point solution, the Pa bosons are fully condensed and 
their fluctuations were ignored for simplicity, the question arises as to how to evaluation the corresponding Green's 
functions in Eq. (|S43|) . Naively, one would simply ignore these fluctuations; but we found that doing so would produce 
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an unphysical result of a nonzero incoherent contribution to the ISF in the nonintcracting hmit. The problem can be 
resolved by realizing that at the saddle point level, the local constraint in Eq. (|S4[) is only satisfied on average, i.e. 
{Qia) = 0. When fluctuations are considered, a consistent condition imposed by the constraint is {TrQ a{T)Q = 
where Qq = (1/^) X^iea '5ia with N the number of a-sublattice sites. Evaluating the latter leads to the following 
expression 



k 



= 0, 



(S48) 



to leading order in the boson correlations. As a result, the Green's function of the boson in AJ^^°'^ given in Eq. ()S43P 
can be expressed in terms of those of the d-e bosons; leading to 

ALr'(T) = E {pL(r.elq(r)g„q(0)) +pL(r.daq(r)dlq(0)) +pt [(Tr4q(r)4q(0)> + {TrL^{T)e^^m] }, (S49) 
q 

where 

Pla=Ploa-Po^ pL=pIo9-Po: ptj = PotPol ' pI', With pi = 2df^ / {pl^ + pl^) . (S50) 

We have verified that the above results hold even in the presence of pairing between the bosons. The normal and 
anomalous Green's function of the fluctuating doublons and holons in Eq. (jS49p are readily evaluated using the boson 
wavefunctions obtained from the Bogoliubov transformation in Eqs. (|S22p . 

s=± 



s=± 
s=± 



(S51) 
(S52) 
(S53) 
(S54) 



s=± 



where the ub^E) = l/[exp (/?£') — 1] is the boson distribution function. Substituting Eqs (|S5HIS54| ) in Eq. (jS49p . 
inserting the result into Eq. (jS45|) . and carrying out the imaginary-time integral, we obtain the final expression for 
the incoherent part of the ISF, 



N. 



r°V) = -E E-'^{ 



k.q.cr s— ± 



PlaKc^s? + pLWccs? + 2ptRe«q,7.„q,)] |i)ak.P<5[c^ -El- (q)]}. (S55) 



Remarkably, at [/ = 0, where poa = c^o = eo = 1/2, p^^^ = Pqo- = P'ta = Pa'a = in Eq. (jS55|) and the incoherent 
spectral function is completely suppressed, recovering the nonintcracting limit. At any finite [/, our numerical calcu- 
lations show that the spectral sum rule is approximately satisfied and an overall normalization factor close to one is 
applied to Eq. ([SSS]) to ensure W'^^""^ + W^°^ = 1 in the plotted ISF figures. The obtained ISF has been shown in 
Fig. 3 of the main text. To illustrate in more detail the transfer of the coherent QP weight to the incoherent part 
of the spectral function with increasing Hubbard U , we plot in Fig. [S4] separately the coherent (shaded red areas), 
incoherent (blue solid lines) , and the total (black solid lines) ISF for several values of the Hubbard U across the Mott 
transition. 
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